Предварительный анализ данных.

1. Описание признаков присутствует в сопроводительном файле. Прочитаем данные. Среди данных не может быть отрицательных значений, NA обозначаются как -999.

library(readxl)
library(dplyr)
library(tidyr)
df <- read_excel("Sleep/SLEEP_shortname.xls")
df[df < 0] <- NA

head(df)
## # A tibble: 6 × 11
##   NAME     BODY_…¹ BRAIN…² SLOWW…³ PARADOX SLEEP LIFES…⁴ GESTT…⁵ PRED_…⁶ EXP_IND
##   <chr>      <dbl>   <dbl>   <dbl>   <dbl> <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
## 1 African… 6654     5712      NA      NA     3.3    38.6     645       3       5
## 2 African…    1        6.6     6.3     2     8.3     4.5      42       3       1
## 3 Arctic …    3.38    44.5    NA      NA    12.5    14        60       1       1
## 4 Arc. gr…    0.92     5.7    NA      NA    16.5    NA        25       5       2
## 5 Asian e… 2547     4603       2.1     1.8   3.9    69       624       3       5
## 6 Baboon     10.6    180.      9.1     0.7   9.8    27       180       4       4
## # … with 1 more variable: DANG_IND <dbl>, and abbreviated variable names
## #   ¹​BODY_WEI, ²​BRAIN_WE, ³​SLOWWAVE, ⁴​LIFESPAN, ⁵​GESTTIME, ⁶​PRED_IND
## # ℹ Use `colnames()` to see all variable names

Количество строк в таблице.

nrow(df)
## [1] 62

2. Признаков немного, поэтому рассматривать будем все.

3. Кроме индексов-факторов, все признаки количественные. Индексы - порядковые признаки. Все количественные признаки непрерывные, но можно заметить дискретизацию при округлении. Проверим это, посмотрев на частоты мод и повторов.

mode_rate <- function(x) {
  x <- x[!is.na(x)]
  u <- unique(x)
  tab <- tabulate(match(x, u))
  max(tab) / length(x)
}

repeat_rate <- function(x) {
  x <- x[!is.na(x)]
  u <- unique(x)
  (length(x) - length(u)) / length(x)
}

Отношение частоты моды к числу элементов:

df %>% summarise(BODY_WEI = round(mode_rate(BODY_WEI), 3), BRAIN_WE = round(mode_rate(BRAIN_WE), 3), 
                SLOWWAVE = round(mode_rate(SLOWWAVE), 3), PARADOX = round(mode_rate(PARADOX), 3), 
                SLEEP = round(mode_rate(SLEEP), 3), LIFESPAN = round(mode_rate(LIFESPAN), 3), 
                GESTTIME = round(mode_rate(GESTTIME), 3))
## # A tibble: 1 × 7
##   BODY_WEI BRAIN_WE SLOWWAVE PARADOX SLEEP LIFESPAN GESTTIME
##      <dbl>    <dbl>    <dbl>   <dbl> <dbl>    <dbl>    <dbl>
## 1    0.032    0.032    0.062    0.06 0.052    0.052    0.052

Отношение частоты повторных элементов к числу элементов:

df %>% summarise(BODY_WEI = round(repeat_rate(BODY_WEI), 3), BRAIN_WE = round(repeat_rate(BRAIN_WE), 3), 
                SLOWWAVE = round(repeat_rate(SLOWWAVE), 3), PARADOX = round(repeat_rate(PARADOX), 3), 
                SLEEP = round(repeat_rate(SLEEP), 3), LIFESPAN = round(repeat_rate(LIFESPAN), 3), 
                GESTTIME = round(repeat_rate(GESTTIME), 3))
## # A tibble: 1 × 7
##   BODY_WEI BRAIN_WE SLOWWAVE PARADOX SLEEP LIFESPAN GESTTIME
##      <dbl>    <dbl>    <dbl>   <dbl> <dbl>    <dbl>    <dbl>
## 1    0.032    0.048    0.188     0.4 0.241     0.19    0.155

Будем считать данные дискретными если отношение частоты моды к числу элементов больше 10%, а отношение частоты повторных элементов к числу элементов больше 15%. Тогда BODY_WEI, BRAIN_WE будут непрерывными, а BRAIN_WE, PARADOX, SLEEP, LIFESPAN, GESTTIME дискретными.

4. Не актуально для текущих данных.

5. Посмотрим на данные.

library(ggplot2)
library(GGally)
df %>% dplyr::select(-NAME) %>%
  ggpairs(diag=list(continuous = "barDiag"), 
          columns = c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))

6. Преобразуем данные. Переведем вес животного в граммы, продолжительность жизни в дни. Логарифмируем данные, чтобы было легче наблюдать корреляции (выше заметны сильно отличающиеся индивиды, это слоны). Факторизуем индексы.

dfNew <- df %>% mutate(PRED_IND = as.factor(PRED_IND), EXP_IND = as.factor(EXP_IND),
                   DANG_IND = as.factor(DANG_IND), 
                   BODY_WEI = log(BODY_WEI * 1000), BRAIN_WE = log(BRAIN_WE), 
                   LIFESPAN = log(LIFESPAN * 365.25), GESTTIME = log(GESTTIME))

Посмотрим на новые данные. Сгруппируем индивидов по факторам.

No factorization

dfNew %>% dplyr::select(-NAME) %>%
  ggpairs(diag=list(continuous = "barDiag"), 
          columns = c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))

PRED_IND

dfNew %>% dplyr::select(-NAME) %>%
  ggpairs(diag=list(continuous = "barDiag"), aes(colour=PRED_IND), legend=1,
          columns=c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))

EXP_IND

dfNew %>% dplyr::select(-NAME) %>%
  ggpairs(diag=list(continuous = "barDiag"), aes(colour=EXP_IND), legend=1,
          columns=c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))

DANG_IND

dfNew %>% dplyr::select(-NAME) %>%
  ggpairs(diag=list(continuous = "barDiag"), aes(colour=DANG_IND), legend=1,
          columns=c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))

7. Аутлайнеров нет.

8. Неоднородности нет.

9. Не актуально для текущих данных.

10. Используем описательные статистики для распределений признаков в новой выборке.

library(moments)
characteristics <- function(x) {
  c(mean = round(mean(x, na.rm = TRUE), 3),
    median = round(median(x, na.rm = TRUE), 3),
    var = round(var(x, na.rm = TRUE), 3),
    skewness = round(skewness(x, na.rm = TRUE), 3),
    kurtosis = round(kurtosis(x, na.rm = TRUE) - 3, 3))
}

data.frame(lapply(as.list(dfNew[2:8]), characteristics)) %>% t()
##            mean median    var skewness kurtosis
## BODY_WEI  8.245  8.114  9.754    0.149   -0.441
## BRAIN_WE  3.140  2.848  5.985    0.045   -0.540
## SLOWWAVE  8.673  8.350 13.443    0.289   -0.334
## PARADOX   1.972  1.800  2.081    1.409    1.976
## SLEEP    10.533 10.450 21.222    0.196   -0.567
## LIFESPAN  8.494  8.613  0.890   -0.167   -0.841
## GESTTIME  4.444  4.360  1.124    0.049   -1.075

О виде распределений и о сравнении распределений.

Индивидуальное задание 2.

1. Проверим распределения признаков на нормальность.

library(reshape)

Изобразим гистограммы признаков на фоне плотности нормального распределения.

dfMelt <- dfNew %>% 
  dplyr::select(-PRED_IND, -DANG_IND) %>% 
  as.data.frame(dfNew) %>% 
  melt(id.vars = c("NAME", "EXP_IND"))

ggplot(dfMelt, aes(x = value)) +
        geom_histogram(aes(y = ..density..)) +
        facet_wrap(~variable, scales = "free") +
        labs(x = "", y = "") +
        geom_line(aes(y = dnorm(value,
                      mean = tapply(value, variable, mean, na.rm = TRUE)[PANEL],
                      sd = tapply(value, variable, sd, na.rm = TRUE)[PANEL]
                      )), color = "blue")

Изобразим выборочную функцию распредления признаков на фоне функции распределения нормального распределения.

ggplot(dfMelt, aes(x = value)) +
        stat_ecdf() +
        facet_wrap(~variable, scales = "free") +
        labs(x = "", y = "") +
        geom_line(aes(y = pnorm(value,
                      mean = tapply(value, variable, mean, na.rm = TRUE)[PANEL],
                      sd = tapply(value, variable, sd, na.rm = TRUE)[PANEL]
                      )), color = "blue")

Давайте еще нарисуем PP-plot и QQ-plot.

library(qqplotr) 
ggplot(dfMelt, aes(sample = value)) +
        stat_pp_point(size = 1) +
        facet_wrap(~variable, scales = "free") +
        labs(x = "", y = "") +
        stat_pp_line(color = "blue")

ggplot(dfMelt, aes(sample = value)) +
        stat_qq_point(size = 1) +
        facet_wrap(~variable, scales = "free") +
        labs(x = "", y = "") +
        stat_qq_line(color = "blue")

2. Похоже, что нормальность имеет место. Проверим это с помощью критериев.

library(nortest)

У нас достаточно повторений в некоторых выборках, чтобы не использовать критерий Пирсона, но все ранво проведем его, скорее, ради академического интереса. Критерий Лиллиефорса - это модификация критерия Колмогорова-Смирнова для проверки сложных гипотез нормальности данных. Критерий Андерона-Дарлинга - это один из критериев типа w^2. Критерий Шапира-Уилка - примерно квадрат корреляции между x и y в n.p.p.

normality.tests <- function(x) {
  c(Pearson = round(pearson.test(x)$p, 5),
    Lilliefors = round(lillie.test(x)$p, 5),
    Anderson.Darling = round(ad.test(x)$p, 5),
    Shapiro.Wilk = round(shapiro.test(x)$p, 5))
}

data.frame(lapply(as.list(dfNew[2:8]), normality.tests)) %>% t()
##          Pearson Lilliefors Anderson.Darling Shapiro.Wilk
## BODY_WEI 0.26277    0.04414          0.46210      0.72717
## BRAIN_WE 0.71196    0.70193          0.64143      0.69703
## SLOWWAVE 0.86095    0.78346          0.90556      0.77036
## PARADOX  0.13013    0.02883          0.00027      0.00009
## SLEEP    0.07080    0.81866          0.37487      0.13963
## LIFESPAN 0.47137    0.45797          0.20207      0.29791
## GESTTIME 0.51083    0.22157          0.11163      0.12164

Гипотеза о нормальности распределения для быстрого сна отвергается многими критериями с уровнем значисости 0.05. Гипотеза нормльности распределения для веса животного отвергается критерием Лиллиефорса с уровнем значисости 0.05. Остальные критерии гипотезы не отвергли с уровнем значисости 0.05.

Индивидуальное задание 1.

3. Опишем разницу между животными по степени опасности места, где они спят.

BODY_WEI

dfNew %>% ggplot(aes(y=BODY_WEI, colour=EXP_IND)) + geom_boxplot()

BRAIN_WE

dfNew %>% ggplot(aes(y=BRAIN_WE, colour=EXP_IND)) + geom_boxplot()

SLOWWAVE

dfNew %>% ggplot(aes(y=SLOWWAVE, colour=EXP_IND)) + geom_boxplot()

PARADOX

dfNew %>% ggplot(aes(y=PARADOX, colour=EXP_IND)) + geom_boxplot()

SLEEP

dfNew %>% ggplot(aes(y=SLEEP, colour=EXP_IND)) + geom_boxplot()

LIFESPAN

dfNew %>% ggplot(aes(y=LIFESPAN, colour=EXP_IND)) + geom_boxplot()

GESTTIME

dfNew %>% ggplot(aes(y=GESTTIME, colour=EXP_IND)) + geom_boxplot()

Можно заметить, что животные, живущие в более защищенных местах, 1,2) имеют меньший вес тела и мозга; 3,4,5) чаще дольше спят в общем и в разных фазах сна по отдельности; 6,7) имеют менее длительные продолжительность жизни и период вынашивания потомства.

Давайте проверим наблюдения 1,3) c помощью критериев для животных, живущих в ниболее и в наименее защищенных местах.

dfNew1 <- dplyr::filter(dfNew, EXP_IND == 1)
dfNew5 <- dplyr::filter(dfNew, EXP_IND == 5)

nrow(dfNew1)
## [1] 27
nrow(dfNew5)
## [1] 13

Ранее мы узнали, что многие признаки распределены нормально, но сейчас мы будем сравнивать разные группы индивидов, признаки внутри которых могут уже не имеют нормальность.

ggplot(dfNew1, aes(sample = BODY_WEI)) +
        stat_qq_point(size = 1) +
        labs(x = "", y = "") +
        stat_qq_line(color = "blue")

ggplot(dfNew1, aes(sample = SLEEP)) +
        stat_qq_point(size = 1) +
        labs(x = "", y = "") +
        stat_qq_line(color = "blue")

ggplot(dfNew5, aes(sample = BODY_WEI)) +
        stat_qq_point(size = 1) +
        labs(x = "", y = "") +
        stat_qq_line(color = "blue")

ggplot(dfNew5, aes(sample = SLEEP)) +
        stat_qq_point(size = 1) +
        labs(x = "", y = "") +
        stat_qq_line(color = "blue")

Используем те же критерии, что использовали ранее.

Вес тела для животных, спящих в более защищенном месте.

normality.tests(dfNew1$BODY_WEI)
##          Pearson       Lilliefors Anderson.Darling     Shapiro.Wilk 
##          0.07419          0.42187          0.30137          0.24588

Сон для животных, спящих в более защищенном месте.

normality.tests(dfNew1$SLEEP)
##          Pearson       Lilliefors Anderson.Darling     Shapiro.Wilk 
##          0.16492          0.16347          0.09438          0.12591

Вес тела для животных, спящих в менее защищенном месте.

normality.tests(dfNew5$BODY_WEI)
##          Pearson       Lilliefors Anderson.Darling     Shapiro.Wilk 
##          0.36851          0.78725          0.91119          0.99223

Сон для животных, спящих в менее защищенном месте.

normality.tests(dfNew5$SLEEP)
##          Pearson       Lilliefors Anderson.Darling     Shapiro.Wilk 
##          0.00889          0.00047          0.00476          0.00526

Нормальность отверглась с уровням значимости 0.05 толко для сная животных, спящих в менее защищенном месте. Причем по всем критериям отверглась.

Итак. Группы у нас независимые и некоторые из них имеют нормальное распределение.

4. Сравним распределения с помощью t-критериев для независимых выборок. Но сначала проверим гипотезу о равенстве дисперсий распределений.

Для веса тела можно использовать критерий Фишера, а для сна нельзя, так как для сна не выполняется нормальность, поэтому ему доверять не будем, но проведем ради академического интереса. Так же используем критерий Левена с модулями.

var.tests <- function(x, y) {
  c(Fisher = round(var.test(x, y)$p.value, 5),
    Levene = round(t.test(abs(x - mean(x, na.rm = TRUE)), 
                          abs(y - mean(y, na.rm = TRUE)))$p.value, 5))
}

Вес тела:

var.tests(dfNew1$BODY_WEI, dfNew5$BODY_WEI)
##  Fisher  Levene 
## 0.44453 0.28662

Сон:

var.tests(dfNew1$SLEEP, dfNew5$SLEEP)
##  Fisher  Levene 
## 0.01559 0.00497

Гипотеза о равных дисперсиях не отверглась для веса, а для сна отверглась с уровнем значимости 0.05.

Посмотрим на t-критерии с равными и неравными дисперсиями для независимых выборок.

t.tests <- function(x, y) {
  c(Two.sample.t.test = t.test(x, y, var.equal=TRUE)$p.value,
    Welch.two.sample.t.test = t.test(x, y, var.equal=FALSE)$p.value)
}

Вес тела:

t.tests(dfNew1$BODY_WEI, dfNew5$BODY_WEI)
##       Two.sample.t.test Welch.two.sample.t.test 
##            6.044937e-07            4.735032e-07

Сон:

t.tests(dfNew1$SLEEP, dfNew5$SLEEP)
##       Two.sample.t.test Welch.two.sample.t.test 
##            1.292774e-07            1.676559e-10

Первый критерий точный, второй используется не по назначению (так как дисперсии разные), третий и четвертный ассимптотические. Как и предполагалось отверглись гипотезы в пользу того, что мы наблюдали на box-plot-ах выше.

5,6. Так как в случае сна удалось проверить гипотезу только асимптотическим критерием, а выборки у нас небольшие, воспользуемся другими критериями. Выбросов у нас нет, но все же используем критерий Уолкоксона ради академического интереса.

Вес:

wilcox.test(dfNew1$BODY_WEI, dfNew5$BODY_WEI, paired=FALSE)$p.value
## [1] 2.183504e-05

Сон:

wilcox.test(dfNew1$SLEEP, dfNew5$SLEEP, paired=FALSE)$p.value
## [1] 1.002996e-05

Мощность у этого критерия меньшая, но гипотезы так же отверглись.

7. Наконец используем двухвыборочный критерий Колмогорова-Смирнова, который имеет меньшую мощность, но у которого более общая альтернативная гипотеза.

Вес:

ks.test(dfNew1$BODY_WEI, dfNew5$BODY_WEI)$p.value
## [1] 0.0001415934

Сон:

ks.test(dfNew1$SLEEP, dfNew5$SLEEP)$p.value
## [1] 4.413664e-05

Гипотезы так же отверглись.

8. Не актуально для текущих данных.

Об анализе зависимостей.

Индивидуальное задание 3.

1. Посмотрим еще раз на ggpairs plot.

No factorization

dfNew %>% dplyr::select(-NAME) %>%
  ggpairs(diag=list(continuous = "barDiag"), 
          columns = c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))

PRED_IND

dfNew %>% dplyr::select(-NAME) %>%
  ggpairs(diag=list(continuous = "barDiag"), aes(colour=PRED_IND), legend=1,
          columns=c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))

EXP_IND

dfNew %>% dplyr::select(-NAME) %>%
  ggpairs(diag=list(continuous = "barDiag"), aes(colour=EXP_IND), legend=1,
          columns=c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))

DANG_IND

dfNew %>% dplyr::select(-NAME) %>%
  ggpairs(diag=list(continuous = "barDiag"), aes(colour=DANG_IND), legend=1,
          columns=c("BODY_WEI", "BRAIN_WE", "SLOWWAVE", "PARADOX", "SLEEP", "LIFESPAN", "GESTTIME"))

Можно заметить 1) положительную корреляцию веса тела, веса мозга, продолжительности жизни и периода вынашивания потомства; 2) положительную корреляцию фаз сна и сна в общем смысле; 3) отрицательную корреляцию всего, что было перечислено в пункте 1), со всем, что было перечислено в пункте 2).

2. Посмотрим на корреляцию Пирсона между признаками.

dplyr::select(dfNew, -NAME) %>%
mutate(PRED_IND = as.numeric(PRED_IND), EXP_IND = as.numeric(EXP_IND), DANG_IND = as.numeric(DANG_IND)) %>%
cor(method = "pearson", use = "pairwise.complete.obs") %>%
melt() %>%
ggplot(aes(X1, X2)) +
  geom_raster(aes(fill = value)) +
  geom_text(aes(label = round(value, 3))) +
  scale_fill_gradient2(low=colors()[555], mid=colors()[1], high=colors()[26]) + 
  ggtitle("dfNew pearson") +
  theme(axis.text.x = element_text(angle = 50, hjust = 1))

Видим подтверждение того, что наблюдали на ggpairs, еще заметны отрицательные корреляции индексов и сна. То есть чем в большей опасности находится животное, тем меньше оно спит. Еще заметна сильная корреляция между размерами животного и незащищенностью места, где оно спит.

3. Теперь посмотрим на корреляции Спирмана между признаками.

dplyr::select(dfNew, -NAME) %>%
mutate(PRED_IND = as.numeric(PRED_IND), EXP_IND = as.numeric(EXP_IND), DANG_IND = as.numeric(DANG_IND)) %>%
cor(method = "spearman", use = "pairwise.complete.obs") %>%
melt() %>%
ggplot(aes(X1, X2)) +
  geom_raster(aes(fill = value)) +
  geom_text(aes(label = round(value, 3))) +
  scale_fill_gradient2(low=colors()[555], mid=colors()[1], high=colors()[26]) + 
  ggtitle("dfNew spearman") +
  theme(axis.text.x = element_text(angle = 50, hjust = 1))

Больших изменений нет (так как у нас нет выбросов и распределения близкие к нормальным).

4. Также посчитаем интересные частные корреляции.

library(ppcor)

Посмотрим на частную корреляцию веса животного и сна за вычетом идексов опасности и на частную корреляцию степени защищенности животного во время сна и сна за вычетом веса и связанных с ним критериев.

Это будет иметь смысл, так как между индексами и весом/сном коэффициент Пирсона значим, то есть есть линейная зависимость между индексами и весом/сном, которая будет вычитаться при подсчете частных корреляций.

Частая корреляция веса животного и сна за вычетом идексов опасности.

((dplyr::select(dfNew, -NAME, -BRAIN_WE, -LIFESPAN, -GESTTIME, -SLOWWAVE, -PARADOX) %>%
  mutate(PRED_IND = as.numeric(PRED_IND), EXP_IND = as.numeric(EXP_IND), DANG_IND = as.numeric(DANG_IND)) %>%
  drop_na() %>%
  pcor(method = "spearman"))$estimate %>%
as.data.frame())["SLEEP", "BODY_WEI"]
## [1] -0.2668484

Отрицательная корреляции между сном и весом животного остается. Животные с большим весом спят меньше.

Часнтая корреляция степени защищенности животного во время сна и сна за вычетом веса и связанных с ним критериев.

((dplyr::select(dfNew, -NAME, -PRED_IND, -DANG_IND, -SLOWWAVE, -PARADOX) %>%
  mutate(EXP_IND = as.numeric(EXP_IND)) %>%
  drop_na() %>%
  pcor(method = "spearman"))$estimate %>% 
as.data.frame())["SLEEP", "EXP_IND"]
## [1] -0.4055336

Отрицательная корреляция между сном и индексом защищенности животного во время сна остается. Животные, находящиеся в большей опасности во время сна, спят меньше.